home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / gnu / calc202a.lha / calc-2.02a / calc-poly.el < prev    next >
Lisp/Scheme  |  1993-06-01  |  36KB  |  1,196 lines
  1. ;; Calculator for GNU Emacs, part II [calc-poly.el]
  2. ;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
  3. ;; Written by Dave Gillespie, daveg@synaptics.com.
  4.  
  5. ;; This file is part of GNU Emacs.
  6.  
  7. ;; GNU Emacs is distributed in the hope that it will be useful,
  8. ;; but WITHOUT ANY WARRANTY.  No author or distributor
  9. ;; accepts responsibility to anyone for the consequences of using it
  10. ;; or for whether it serves any particular purpose or works at all,
  11. ;; unless he says so in writing.  Refer to the GNU Emacs General Public
  12. ;; License for full details.
  13.  
  14. ;; Everyone is granted permission to copy, modify and redistribute
  15. ;; GNU Emacs, but only under the conditions described in the
  16. ;; GNU Emacs General Public License.   A copy of this license is
  17. ;; supposed to have been given to you along with GNU Emacs so you
  18. ;; can know your rights and responsibilities.  It should be in a
  19. ;; file named COPYING.  Among other things, the copyright notice
  20. ;; and this notice must be preserved on all copies.
  21.  
  22.  
  23.  
  24. ;; This file is autoloaded from calc-ext.el.
  25. (require 'calc-ext)
  26.  
  27. (require 'calc-macs)
  28.  
  29. (defun calc-Need-calc-poly () nil)
  30.  
  31.  
  32. (defun calcFunc-pcont (expr &optional var)
  33.   (cond ((Math-primp expr)
  34.      (cond ((Math-zerop expr) 1)
  35.            ((Math-messy-integerp expr) (math-trunc expr))
  36.            ((Math-objectp expr) expr)
  37.            ((or (equal expr var) (not var)) 1)
  38.            (t expr)))
  39.     ((eq (car expr) '*)
  40.      (math-mul (calcFunc-pcont (nth 1 expr) var)
  41.            (calcFunc-pcont (nth 2 expr) var)))
  42.     ((eq (car expr) '/)
  43.      (math-div (calcFunc-pcont (nth 1 expr) var)
  44.            (calcFunc-pcont (nth 2 expr) var)))
  45.     ((and (eq (car expr) '^) (Math-natnump (nth 2 expr)))
  46.      (math-pow (calcFunc-pcont (nth 1 expr) var) (nth 2 expr)))
  47.     ((memq (car expr) '(neg polar))
  48.      (calcFunc-pcont (nth 1 expr) var))
  49.     ((consp var)
  50.      (let ((p (math-is-polynomial expr var)))
  51.        (if p
  52.            (let ((lead (nth (1- (length p)) p))
  53.              (cont (math-poly-gcd-list p)))
  54.          (if (math-guess-if-neg lead)
  55.              (math-neg cont)
  56.            cont))
  57.          1)))
  58.     ((memq (car expr) '(+ - cplx sdev))
  59.      (let ((cont (calcFunc-pcont (nth 1 expr) var)))
  60.        (if (eq cont 1)
  61.            1
  62.          (let ((c2 (calcFunc-pcont (nth 2 expr) var)))
  63.            (if (and (math-negp cont)
  64.             (if (eq (car expr) '-) (math-posp c2) (math-negp c2)))
  65.            (math-neg (math-poly-gcd cont c2))
  66.          (math-poly-gcd cont c2))))))
  67.     (var expr)
  68.     (t 1))
  69. )
  70.  
  71. (defun calcFunc-pprim (expr &optional var)
  72.   (let ((cont (calcFunc-pcont expr var)))
  73.     (if (math-equal-int cont 1)
  74.     expr
  75.       (math-poly-div-exact expr cont var)))
  76. )
  77.  
  78. (defun math-div-poly-const (expr c)
  79.   (cond ((memq (car-safe expr) '(+ -))
  80.      (list (car expr)
  81.            (math-div-poly-const (nth 1 expr) c)
  82.            (math-div-poly-const (nth 2 expr) c)))
  83.     (t (math-div expr c)))
  84. )
  85.  
  86. (defun calcFunc-pdeg (expr &optional var)
  87.   (if (Math-zerop expr)
  88.       '(neg (var inf var-inf))
  89.     (if var
  90.     (or (math-polynomial-p expr var)
  91.         (math-reject-arg expr "Expected a polynomial"))
  92.       (math-poly-degree expr)))
  93. )
  94.  
  95. (defun math-poly-degree (expr)
  96.   (cond ((Math-primp expr)
  97.      (if (eq (car-safe expr) 'var) 1 0))
  98.     ((eq (car expr) 'neg)
  99.      (math-poly-degree (nth 1 expr)))
  100.     ((eq (car expr) '*)
  101.      (+ (math-poly-degree (nth 1 expr))
  102.         (math-poly-degree (nth 2 expr))))
  103.     ((eq (car expr) '/)
  104.      (- (math-poly-degree (nth 1 expr))
  105.         (math-poly-degree (nth 2 expr))))
  106.     ((and (eq (car expr) '^) (natnump (nth 2 expr)))
  107.      (* (math-poly-degree (nth 1 expr)) (nth 2 expr)))
  108.     ((memq (car expr) '(+ -))
  109.      (max (math-poly-degree (nth 1 expr))
  110.           (math-poly-degree (nth 2 expr))))
  111.     (t 1))
  112. )
  113.  
  114. (defun calcFunc-plead (expr var)
  115.   (cond ((eq (car-safe expr) '*)
  116.      (math-mul (calcFunc-plead (nth 1 expr) var)
  117.            (calcFunc-plead (nth 2 expr) var)))
  118.     ((eq (car-safe expr) '/)
  119.      (math-div (calcFunc-plead (nth 1 expr) var)
  120.            (calcFunc-plead (nth 2 expr) var)))
  121.     ((and (eq (car-safe expr) '^) (math-natnump (nth 2 expr)))
  122.      (math-pow (calcFunc-plead (nth 1 expr) var) (nth 2 expr)))
  123.     ((Math-primp expr)
  124.      (if (equal expr var)
  125.          1
  126.        expr))
  127.     (t
  128.      (let ((p (math-is-polynomial expr var)))
  129.        (if (cdr p)
  130.            (nth (1- (length p)) p)
  131.          1))))
  132. )
  133.  
  134.  
  135.  
  136.  
  137.  
  138. ;;; Polynomial quotient, remainder, and GCD.
  139. ;;; Originally by Ove Ewerlid (ewerlid@mizar.DoCS.UU.SE).
  140. ;;; Modifications and simplifications by daveg.
  141.  
  142. (setq math-poly-modulus 1)
  143.  
  144. ;;; Return gcd of two polynomials
  145. (defun calcFunc-pgcd (pn pd)
  146.   (if (math-any-floats pn)
  147.       (math-reject-arg pn "Coefficients must be rational"))
  148.   (if (math-any-floats pd)
  149.       (math-reject-arg pd "Coefficients must be rational"))
  150.   (let ((calc-prefer-frac t)
  151.     (math-poly-modulus (math-poly-modulus pn pd)))
  152.     (math-poly-gcd pn pd))
  153. )
  154.  
  155. ;;; Return only quotient to top of stack (nil if zero)
  156. (defun calcFunc-pdiv (pn pd &optional base)
  157.   (let* ((calc-prefer-frac t)
  158.      (math-poly-modulus (math-poly-modulus pn pd))
  159.      (res (math-poly-div pn pd base)))
  160.     (setq calc-poly-div-remainder (cdr res))
  161.     (car res))
  162. )
  163.  
  164. ;;; Return only remainder to top of stack
  165. (defun calcFunc-prem (pn pd &optional base)
  166.   (let ((calc-prefer-frac t)
  167.     (math-poly-modulus (math-poly-modulus pn pd)))
  168.     (cdr (math-poly-div pn pd base)))
  169. )
  170.  
  171. (defun calcFunc-pdivrem (pn pd &optional base)
  172.   (let* ((calc-prefer-frac t)
  173.      (math-poly-modulus (math-poly-modulus pn pd))
  174.      (res (math-poly-div pn pd base)))
  175.     (list 'vec (car res) (cdr res)))
  176. )
  177.  
  178. (defun calcFunc-pdivide (pn pd &optional base)
  179.   (let* ((calc-prefer-frac t)
  180.      (math-poly-modulus (math-poly-modulus pn pd))
  181.      (res (math-poly-div pn pd base)))
  182.     (math-add (car res) (math-div (cdr res) pd)))
  183. )
  184.  
  185.  
  186. ;;; Multiply two terms, expanding out products of sums.
  187. (defun math-mul-thru (lhs rhs)
  188.   (if (memq (car-safe lhs) '(+ -))
  189.       (list (car lhs)
  190.         (math-mul-thru (nth 1 lhs) rhs)
  191.         (math-mul-thru (nth 2 lhs) rhs))
  192.     (if (memq (car-safe rhs) '(+ -))
  193.     (list (car rhs)
  194.           (math-mul-thru lhs (nth 1 rhs))
  195.           (math-mul-thru lhs (nth 2 rhs)))
  196.       (math-mul lhs rhs)))
  197. )
  198.  
  199. (defun math-div-thru (num den)
  200.   (if (memq (car-safe num) '(+ -))
  201.       (list (car num)
  202.         (math-div-thru (nth 1 num) den)
  203.         (math-div-thru (nth 2 num) den))
  204.     (math-div num den))
  205. )
  206.  
  207.  
  208. ;;; Sort the terms of a sum into canonical order.
  209. (defun math-sort-terms (expr)
  210.   (if (memq (car-safe expr) '(+ -))
  211.       (math-list-to-sum
  212.        (sort (math-sum-to-list expr)
  213.          (function (lambda (a b) (math-beforep (car a) (car b))))))
  214.     expr)
  215. )
  216.  
  217. (defun math-list-to-sum (lst)
  218.   (if (cdr lst)
  219.       (list (if (cdr (car lst)) '- '+)
  220.         (math-list-to-sum (cdr lst))
  221.         (car (car lst)))
  222.     (if (cdr (car lst))
  223.     (math-neg (car (car lst)))
  224.       (car (car lst))))
  225. )
  226.  
  227. (defun math-sum-to-list (tree &optional neg)
  228.   (cond ((eq (car-safe tree) '+)
  229.      (nconc (math-sum-to-list (nth 1 tree) neg)
  230.         (math-sum-to-list (nth 2 tree) neg)))
  231.     ((eq (car-safe tree) '-)
  232.      (nconc (math-sum-to-list (nth 1 tree) neg)
  233.         (math-sum-to-list (nth 2 tree) (not neg))))
  234.     (t (list (cons tree neg))))
  235. )
  236.  
  237. ;;; Check if the polynomial coefficients are modulo forms.
  238. (defun math-poly-modulus (expr &optional expr2)
  239.   (or (math-poly-modulus-rec expr)
  240.       (and expr2 (math-poly-modulus-rec expr2))
  241.       1)
  242. )
  243.  
  244. (defun math-poly-modulus-rec (expr)
  245.   (if (and (eq (car-safe expr) 'mod) (Math-natnump (nth 2 expr)))
  246.       (list 'mod 1 (nth 2 expr))
  247.     (and (memq (car-safe expr) '(+ - * /))
  248.      (or (math-poly-modulus-rec (nth 1 expr))
  249.          (math-poly-modulus-rec (nth 2 expr)))))
  250. )
  251.  
  252.  
  253. ;;; Divide two polynomials.  Return (quotient . remainder).
  254. (defun math-poly-div (u v &optional math-poly-div-base)
  255.   (if math-poly-div-base
  256.       (math-do-poly-div u v)
  257.     (math-do-poly-div (calcFunc-expand u) (calcFunc-expand v)))
  258. )
  259. (setq math-poly-div-base nil)
  260.  
  261. (defun math-poly-div-exact (u v &optional base)
  262.   (let ((res (math-poly-div u v base)))
  263.     (if (eq (cdr res) 0)
  264.     (car res)
  265.       (math-reject-arg (list 'vec u v) "Argument is not a polynomial")))
  266. )
  267.  
  268. (defun math-do-poly-div (u v)
  269.   (cond ((math-constp u)
  270.      (if (math-constp v)
  271.          (cons (math-div u v) 0)
  272.        (cons 0 u)))
  273.     ((math-constp v)
  274.      (cons (if (eq v 1)
  275.            u
  276.          (if (memq (car-safe u) '(+ -))
  277.              (math-add-or-sub (math-poly-div-exact (nth 1 u) v)
  278.                       (math-poly-div-exact (nth 2 u) v)
  279.                       nil (eq (car u) '-))
  280.            (math-div u v)))
  281.            0))
  282.     ((Math-equal u v)
  283.      (cons math-poly-modulus 0))
  284.     ((and (math-atomic-factorp u) (math-atomic-factorp v))
  285.      (cons (math-simplify (math-div u v)) 0))
  286.     (t
  287.      (let ((base (or math-poly-div-base
  288.              (math-poly-div-base u v)))
  289.            vp up res)
  290.        (if (or (null base)
  291.            (null (setq vp (math-is-polynomial v base nil 'gen))))
  292.            (cons 0 u)
  293.          (setq up (math-is-polynomial u base nil 'gen)
  294.            res (math-poly-div-coefs up vp))
  295.          (cons (math-build-polynomial-expr (car res) base)
  296.            (math-build-polynomial-expr (cdr res) base))))))
  297. )
  298.  
  299. (defun math-poly-div-rec (u v)
  300.   (cond ((math-constp u)
  301.      (math-div u v))
  302.     ((math-constp v)
  303.      (if (eq v 1)
  304.          u
  305.        (if (memq (car-safe u) '(+ -))
  306.            (math-add-or-sub (math-poly-div-rec (nth 1 u) v)
  307.                 (math-poly-div-rec (nth 2 u) v)
  308.                 nil (eq (car u) '-))
  309.          (math-div u v))))
  310.     ((Math-equal u v) math-poly-modulus)
  311.     ((and (math-atomic-factorp u) (math-atomic-factorp v))
  312.      (math-simplify (math-div u v)))
  313.     (math-poly-div-base
  314.      (math-div u v))
  315.     (t
  316.      (let ((base (math-poly-div-base u v))
  317.            vp up res)
  318.        (if (or (null base)
  319.            (null (setq vp (math-is-polynomial v base nil 'gen))))
  320.            (math-div u v)
  321.          (setq up (math-is-polynomial u base nil 'gen)
  322.            res (math-poly-div-coefs up vp))
  323.          (math-add (math-build-polynomial-expr (car res) base)
  324.                (math-div (math-build-polynomial-expr (cdr res) base)
  325.                  v))))))
  326. )
  327.  
  328. ;;; Divide two polynomials in coefficient-list form.  Return (quot . rem).
  329. (defun math-poly-div-coefs (u v)
  330.   (cond ((null v) (math-reject-arg nil "Division by zero"))
  331.     ((< (length u) (length v)) (cons nil u))
  332.     ((cdr u)
  333.      (let ((q nil)
  334.            (urev (reverse u))
  335.            (vrev (reverse v)))
  336.        (while
  337.            (let ((qk (math-poly-div-rec (math-simplify (car urev))
  338.                         (car vrev)))
  339.              (up urev)
  340.              (vp vrev))
  341.          (if (or q (not (math-zerop qk)))
  342.              (setq q (cons qk q)))
  343.          (while (setq up (cdr up) vp (cdr vp))
  344.            (setcar up (math-sub (car up) (math-mul-thru qk (car vp)))))
  345.          (setq urev (cdr urev))
  346.          up))
  347.        (while (and urev (Math-zerop (car urev)))
  348.          (setq urev (cdr urev)))
  349.        (cons q (nreverse (mapcar 'math-simplify urev)))))
  350.     (t
  351.      (cons (list (math-poly-div-rec (car u) (car v)))
  352.            nil)))
  353. )
  354.  
  355. ;;; Perform a pseudo-division of polynomials.  (See Knuth section 4.6.1.)
  356. ;;; This returns only the remainder from the pseudo-division.
  357. (defun math-poly-pseudo-div (u v)
  358.   (cond ((null v) nil)
  359.     ((< (length u) (length v)) u)
  360.     ((or (cdr u) (cdr v))
  361.      (let ((urev (reverse u))
  362.            (vrev (reverse v))
  363.            up)
  364.        (while
  365.            (let ((vp vrev))
  366.          (setq up urev)
  367.          (while (setq up (cdr up) vp (cdr vp))
  368.            (setcar up (math-sub (math-mul-thru (car vrev) (car up))
  369.                     (math-mul-thru (car urev) (car vp)))))
  370.          (setq urev (cdr urev))
  371.          up)
  372.          (while up
  373.            (setcar up (math-mul-thru (car vrev) (car up)))
  374.            (setq up (cdr up))))
  375.        (while (and urev (Math-zerop (car urev)))
  376.          (setq urev (cdr urev)))
  377.        (nreverse (mapcar 'math-simplify urev))))
  378.     (t nil))
  379. )
  380.  
  381. ;;; Compute the GCD of two multivariate polynomials.
  382. (defun math-poly-gcd (u v)
  383.   (cond ((Math-equal u v) u)
  384.     ((math-constp u)
  385.      (if (Math-zerop u)
  386.          v
  387.        (calcFunc-gcd u (calcFunc-pcont v))))
  388.     ((math-constp v)
  389.      (if (Math-zerop v)
  390.          v
  391.        (calcFunc-gcd v (calcFunc-pcont u))))
  392.     (t
  393.      (let ((base (math-poly-gcd-base u v)))
  394.        (if base
  395.            (math-simplify
  396.         (calcFunc-expand
  397.          (math-build-polynomial-expr
  398.           (math-poly-gcd-coefs (math-is-polynomial u base nil 'gen)
  399.                        (math-is-polynomial v base nil 'gen))
  400.           base)))
  401.          (calcFunc-gcd (calcFunc-pcont u) (calcFunc-pcont u))))))
  402. )
  403.  
  404. (defun math-poly-div-list (lst a)
  405.   (if (eq a 1)
  406.       lst
  407.     (if (eq a -1)
  408.     (math-mul-list lst a)
  409.       (mapcar (function (lambda (x) (math-poly-div-exact x a))) lst)))
  410. )
  411.  
  412. (defun math-mul-list (lst a)
  413.   (if (eq a 1)
  414.       lst
  415.     (if (eq a -1)
  416.     (mapcar 'math-neg lst)
  417.       (and (not (eq a 0))
  418.        (mapcar (function (lambda (x) (math-mul x a))) lst))))
  419. )
  420.  
  421. ;;; Run GCD on all elements in a list.
  422. (defun math-poly-gcd-list (lst)
  423.   (if (or (memq 1 lst) (memq -1 lst))
  424.       (math-poly-gcd-frac-list lst)
  425.     (let ((gcd (car lst)))
  426.       (while (and (setq lst (cdr lst)) (not (eq gcd 1)))
  427.     (or (eq (car lst) 0)
  428.         (setq gcd (math-poly-gcd gcd (car lst)))))
  429.       (if lst (setq lst (math-poly-gcd-frac-list lst)))
  430.       gcd))
  431. )
  432.  
  433. (defun math-poly-gcd-frac-list (lst)
  434.   (while (and lst (not (eq (car-safe (car lst)) 'frac)))
  435.     (setq lst (cdr lst)))
  436.   (if lst
  437.       (let ((denom (nth 2 (car lst))))
  438.     (while (setq lst (cdr lst))
  439.       (if (eq (car-safe (car lst)) 'frac)
  440.           (setq denom (calcFunc-lcm denom (nth 2 (car lst))))))
  441.     (list 'frac 1 denom))
  442.     1)
  443. )
  444.  
  445. ;;; Compute the GCD of two monovariate polynomial lists.
  446. ;;; Knuth section 4.6.1, algorithm C.
  447. (defun math-poly-gcd-coefs (u v)
  448.   (let ((d (math-poly-gcd (math-poly-gcd-list u)
  449.               (math-poly-gcd-list v)))
  450.     (g 1) (h 1) (z 0) hh r delta ghd)
  451.     (while (and u v (Math-zerop (car u)) (Math-zerop (car v)))
  452.       (setq u (cdr u) v (cdr v) z (1+ z)))
  453.     (or (eq d 1)
  454.     (setq u (math-poly-div-list u d)
  455.           v (math-poly-div-list v d)))
  456.     (while (progn
  457.          (setq delta (- (length u) (length v)))
  458.          (if (< delta 0)
  459.          (setq r u u v v r delta (- delta)))
  460.          (setq r (math-poly-pseudo-div u v))
  461.          (cdr r))
  462.       (setq u v
  463.         v (math-poly-div-list r (math-mul g (math-pow h delta)))
  464.         g (nth (1- (length u)) u)
  465.         h (if (<= delta 1)
  466.           (math-mul (math-pow g delta) (math-pow h (- 1 delta)))
  467.         (math-poly-div-exact (math-pow g delta)
  468.                      (math-pow h (1- delta))))))
  469.     (setq v (if r
  470.         (list d)
  471.           (math-mul-list (math-poly-div-list v (math-poly-gcd-list v)) d)))
  472.     (if (math-guess-if-neg (nth (1- (length v)) v))
  473.     (setq v (math-mul-list v -1)))
  474.     (while (>= (setq z (1- z)) 0)
  475.       (setq v (cons 0 v)))
  476.     v)
  477. )
  478.  
  479.  
  480. ;;; Return true if is a factor containing no sums or quotients.
  481. (defun math-atomic-factorp (expr)
  482.   (cond ((eq (car-safe expr) '*)
  483.      (and (math-atomic-factorp (nth 1 expr))
  484.           (math-atomic-factorp (nth 2 expr))))
  485.     ((memq (car-safe expr) '(+ - /))
  486.      nil)
  487.     ((memq (car-safe expr) '(^ neg))
  488.      (math-atomic-factorp (nth 1 expr)))
  489.     (t t))
  490. )
  491.  
  492. ;;; Find a suitable base for dividing a by b.
  493. ;;; The base must exist in both expressions.
  494. ;;; The degree in the numerator must be higher or equal than the
  495. ;;; degree in the denominator.
  496. ;;; If the above conditions are not met the quotient is just a remainder.
  497. ;;; Return nil if this is the case.
  498.  
  499. (defun math-poly-div-base (a b)
  500.   (let (a-base b-base)
  501.     (and (setq a-base (math-total-polynomial-base a))
  502.      (setq b-base (math-total-polynomial-base b))
  503.      (catch 'return
  504.        (while a-base
  505.          (let ((maybe (assoc (car (car a-base)) b-base)))
  506.            (if maybe
  507.            (if (>= (nth 1 (car a-base)) (nth 1 maybe))
  508.                (throw 'return (car (car a-base))))))
  509.          (setq a-base (cdr a-base))))))
  510. )
  511.  
  512. ;;; Same as above but for gcd algorithm.
  513. ;;; Here there is no requirement that degree(a) > degree(b).
  514. ;;; Take the base that has the highest degree considering both a and b.
  515. ;;; ("a^20+b^21+x^3+a+b", "a+b^2+x^5+a^22+b^10") --> (a 22)
  516.  
  517. (defun math-poly-gcd-base (a b)
  518.   (let (a-base b-base)
  519.     (and (setq a-base (math-total-polynomial-base a))
  520.      (setq b-base (math-total-polynomial-base b))
  521.      (catch 'return
  522.        (while (and a-base b-base)
  523.          (if (> (nth 1 (car a-base)) (nth 1 (car b-base)))
  524.          (if (assoc (car (car a-base)) b-base)
  525.              (throw 'return (car (car a-base)))
  526.            (setq a-base (cdr a-base)))
  527.            (if (assoc (car (car b-base)) a-base)
  528.            (throw 'return (car (car b-base)))
  529.          (setq b-base (cdr b-base))))))))
  530. )
  531.  
  532. ;;; Sort a list of polynomial bases.
  533. (defun math-sort-poly-base-list (lst)
  534.   (sort lst (function (lambda (a b)
  535.             (or (> (nth 1 a) (nth 1 b))
  536.                 (and (= (nth 1 a) (nth 1 b))
  537.                  (math-beforep (car a) (car b)))))))
  538. )
  539.  
  540. ;;; Given an expression find all variables that are polynomial bases.
  541. ;;; Return list in the form '( (var1 degree1) (var2 degree2) ... ).
  542. ;;; Note dynamic scope of mpb-total-base.
  543. (defun math-total-polynomial-base (expr)
  544.   (let ((mpb-total-base nil))
  545.     (math-polynomial-base expr 'math-polynomial-p1)
  546.     (math-sort-poly-base-list mpb-total-base))
  547. )
  548.  
  549. (defun math-polynomial-p1 (subexpr)
  550.   (or (assoc subexpr mpb-total-base)
  551.       (memq (car subexpr) '(+ - * / neg))
  552.       (and (eq (car subexpr) '^) (natnump (nth 2 subexpr)))
  553.       (let* ((math-poly-base-variable subexpr)
  554.          (exponent (math-polynomial-p mpb-top-expr subexpr)))
  555.     (if exponent
  556.         (setq mpb-total-base (cons (list subexpr exponent)
  557.                        mpb-total-base)))))
  558.   nil
  559. )
  560.  
  561.  
  562.  
  563.  
  564. (defun calcFunc-factors (expr &optional var)
  565.   (let ((math-factored-vars (if var t nil))
  566.     (math-to-list t)
  567.     (calc-prefer-frac t))
  568.     (or var
  569.     (setq var (math-polynomial-base expr)))
  570.     (let ((res (math-factor-finish
  571.         (or (catch 'factor (math-factor-expr-try var))
  572.             expr))))
  573.       (math-simplify (if (math-vectorp res)
  574.              res
  575.                (list 'vec (list 'vec res 1))))))
  576. )
  577.  
  578. (defun calcFunc-factor (expr &optional var)
  579.   (let ((math-factored-vars nil)
  580.     (math-to-list nil)
  581.     (calc-prefer-frac t))
  582.     (math-simplify (math-factor-finish
  583.             (if var
  584.             (let ((math-factored-vars t))
  585.               (or (catch 'factor (math-factor-expr-try var)) expr))
  586.               (math-factor-expr expr)))))
  587. )
  588.  
  589. (defun math-factor-finish (x)
  590.   (if (Math-primp x)
  591.       x
  592.     (if (eq (car x) 'calcFunc-Fac-Prot)
  593.     (math-factor-finish (nth 1 x))
  594.       (cons (car x) (mapcar 'math-factor-finish (cdr x)))))
  595. )
  596.  
  597. (defun math-factor-protect (x)
  598.   (if (memq (car-safe x) '(+ -))
  599.       (list 'calcFunc-Fac-Prot x)
  600.     x)
  601. )
  602.  
  603. (defun math-factor-expr (expr)
  604.   (cond ((eq math-factored-vars t) expr)
  605.     ((or (memq (car-safe expr) '(* / ^ neg))
  606.          (assq (car-safe expr) calc-tweak-eqn-table))
  607.      (cons (car expr) (mapcar 'math-factor-expr (cdr expr))))
  608.     ((memq (car-safe expr) '(+ -))
  609.      (let* ((math-factored-vars math-factored-vars)
  610.         (y (catch 'factor (math-factor-expr-part expr))))
  611.        (if y
  612.            (math-factor-expr y)
  613.          expr)))
  614.     (t expr))
  615. )
  616.  
  617. (defun math-factor-expr-part (x)    ; uses "expr"
  618.   (if (memq (car-safe x) '(+ - * / ^ neg))
  619.       (while (setq x (cdr x))
  620.     (math-factor-expr-part (car x)))
  621.     (and (not (Math-objvecp x))
  622.      (not (assoc x math-factored-vars))
  623.      (> (math-factor-contains expr x) 1)
  624.      (setq math-factored-vars (cons (list x) math-factored-vars))
  625.      (math-factor-expr-try x)))
  626. )
  627.  
  628. (defun math-factor-expr-try (x)
  629.   (if (eq (car-safe expr) '*)
  630.       (let ((res1 (catch 'factor (let ((expr (nth 1 expr)))
  631.                    (math-factor-expr-try x))))
  632.         (res2 (catch 'factor (let ((expr (nth 2 expr)))
  633.                    (math-factor-expr-try x)))))
  634.     (and (or res1 res2)
  635.          (throw 'factor (math-accum-factors (or res1 (nth 1 expr)) 1
  636.                         (or res2 (nth 2 expr))))))
  637.     (let* ((p (math-is-polynomial expr x 30 'gen))
  638.        (math-poly-modulus (math-poly-modulus expr))
  639.        res)
  640.       (and (cdr p)
  641.        (setq res (math-factor-poly-coefs p))
  642.        (throw 'factor res))))
  643. )
  644.  
  645. (defun math-accum-factors (fac pow facs)
  646.   (if math-to-list
  647.       (if (math-vectorp fac)
  648.       (progn
  649.         (while (setq fac (cdr fac))
  650.           (setq facs (math-accum-factors (nth 1 (car fac))
  651.                          (* pow (nth 2 (car fac)))
  652.                          facs)))
  653.         facs)
  654.     (if (and (eq (car-safe fac) '^) (natnump (nth 2 fac)))
  655.         (setq pow (* pow (nth 2 fac))
  656.           fac (nth 1 fac)))
  657.     (if (eq fac 1)
  658.         facs
  659.       (or (math-vectorp facs)
  660.           (setq facs (if (eq facs 1) '(vec)
  661.                (list 'vec (list 'vec facs 1)))))
  662.       (let ((found facs))
  663.         (while (and (setq found (cdr found))
  664.             (not (equal fac (nth 1 (car found))))))
  665.         (if found
  666.         (progn
  667.           (setcar (cdr (cdr (car found))) (+ pow (nth 2 (car found))))
  668.           facs)
  669.           ;; Put constant term first.
  670.           (if (and (cdr facs) (Math-ratp (nth 1 (nth 1 facs))))
  671.           (cons 'vec (cons (nth 1 facs) (cons (list 'vec fac pow)
  672.                               (cdr (cdr facs)))))
  673.         (cons 'vec (cons (list 'vec fac pow) (cdr facs))))))))
  674.     (math-mul (math-pow fac pow) facs))
  675. )
  676.  
  677. (defun math-factor-poly-coefs (p &optional square-free)    ; uses "x"
  678.   (let (t1 t2)
  679.     (cond ((not (cdr p))
  680.        (or (car p) 0))
  681.  
  682.       ;; Strip off multiples of x.
  683.       ((Math-zerop (car p))
  684.        (let ((z 0))
  685.          (while (and p (Math-zerop (car p)))
  686.            (setq z (1+ z) p (cdr p)))
  687.          (if (cdr p)
  688.          (setq p (math-factor-poly-coefs p square-free))
  689.            (setq p (math-sort-terms (math-factor-expr (car p)))))
  690.          (math-accum-factors x z (math-factor-protect p))))
  691.  
  692.       ;; Factor out content.
  693.       ((and (not square-free)
  694.         (not (eq 1 (setq t1 (math-mul (math-poly-gcd-list p)
  695.                           (if (math-guess-if-neg
  696.                            (nth (1- (length p)) p))
  697.                           -1 1))))))
  698.        (math-accum-factors t1 1 (math-factor-poly-coefs
  699.                      (math-poly-div-list p t1) 'cont)))
  700.  
  701.       ;; Check if linear in x.
  702.       ((not (cdr (cdr p)))
  703.        (math-add (math-factor-protect
  704.               (math-sort-terms
  705.                (math-factor-expr (car p))))
  706.              (math-mul x (math-factor-protect
  707.                   (math-sort-terms
  708.                    (math-factor-expr (nth 1 p)))))))
  709.  
  710.       ;; If symbolic coefficients, use FactorRules.
  711.       ((let ((pp p))
  712.          (while (and pp (or (Math-ratp (car pp))
  713.                 (and (eq (car (car pp)) 'mod)
  714.                      (Math-integerp (nth 1 (car pp)))
  715.                      (Math-integerp (nth 2 (car pp))))))
  716.            (setq pp (cdr pp)))
  717.          pp)
  718.        (let ((res (math-rewrite
  719.                (list 'calcFunc-thecoefs x (cons 'vec p))
  720.                '(var FactorRules var-FactorRules))))
  721.          (or (and (eq (car-safe res) 'calcFunc-thefactors)
  722.               (= (length res) 3)
  723.               (math-vectorp (nth 2 res))
  724.               (let ((facs 1)
  725.                 (vec (nth 2 res)))
  726.             (while (setq vec (cdr vec))
  727.               (setq facs (math-accum-factors (car vec) 1 facs)))
  728.             facs))
  729.          (math-build-polynomial-expr p x))))
  730.  
  731.       ;; Check if rational coefficients (i.e., not modulo a prime).
  732.       ((eq math-poly-modulus 1)
  733.  
  734.        ;; Check if there are any squared terms, or a content not = 1.
  735.        (if (or (eq square-free t)
  736.            (equal (setq t1 (math-poly-gcd-coefs
  737.                     p (setq t2 (math-poly-deriv-coefs p))))
  738.               '(1)))
  739.  
  740.            ;; We now have a square-free polynomial with integer coefs.
  741.            ;; For now, we use a kludgey method that finds linear and
  742.            ;; quadratic terms using floating-point root-finding.
  743.            (if (setq t1 (let ((calc-symbolic-mode nil))
  744.                   (math-poly-all-roots nil p t)))
  745.            (let ((roots (car t1))
  746.              (csign (if (math-negp (nth (1- (length p)) p)) -1 1))
  747.              (expr 1)
  748.              (unfac (nth 1 t1))
  749.              (scale (nth 2 t1)))
  750.              (while roots
  751.                (let ((coef0 (car (car roots)))
  752.                  (coef1 (cdr (car roots))))
  753.              (setq expr (math-accum-factors
  754.                      (if coef1
  755.                      (let ((den (math-lcm-denoms
  756.                              coef0 coef1)))
  757.                        (setq scale (math-div scale den))
  758.                        (math-add
  759.                         (math-add
  760.                          (math-mul den (math-pow x 2))
  761.                          (math-mul (math-mul coef1 den) x))
  762.                         (math-mul coef0 den)))
  763.                        (let ((den (math-lcm-denoms coef0)))
  764.                      (setq scale (math-div scale den))
  765.                      (math-add (math-mul den x)
  766.                            (math-mul coef0 den))))
  767.                      1 expr)
  768.                    roots (cdr roots))))
  769.              (setq expr (math-accum-factors
  770.                  expr 1
  771.                  (math-mul csign
  772.                        (math-build-polynomial-expr
  773.                         (math-mul-list (nth 1 t1) scale)
  774.                         x)))))
  775.          (math-build-polynomial-expr p x))   ; can't factor it.
  776.  
  777.          ;; Separate out the squared terms (Knuth exercise 4.6.2-34).
  778.          ;; This step also divides out the content of the polynomial.
  779.          (let* ((cabs (math-poly-gcd-list p))
  780.             (csign (if (math-negp (nth (1- (length p)) p)) -1 1))
  781.             (t1s (math-mul-list t1 csign))
  782.             (uu nil)
  783.             (v (car (math-poly-div-coefs p t1s)))
  784.             (w (car (math-poly-div-coefs t2 t1s))))
  785.            (while
  786.            (not (math-poly-zerop
  787.              (setq t2 (math-poly-simplify
  788.                    (math-poly-mix
  789.                     w 1 (math-poly-deriv-coefs v) -1)))))
  790.          (setq t1 (math-poly-gcd-coefs v t2)
  791.                uu (cons t1 uu)
  792.                v (car (math-poly-div-coefs v t1))
  793.                w (car (math-poly-div-coefs t2 t1))))
  794.            (setq t1 (length uu)
  795.              t2 (math-accum-factors (math-factor-poly-coefs v t)
  796.                         (1+ t1) 1))
  797.            (while uu
  798.          (setq t2 (math-accum-factors (math-factor-poly-coefs
  799.                            (car uu) t)
  800.                           t1 t2)
  801.                t1 (1- t1)
  802.                uu (cdr uu)))
  803.            (math-accum-factors (math-mul cabs csign) 1 t2))))
  804.  
  805.       ;; Factoring modulo a prime.
  806.       ((and (= (length (setq temp (math-poly-gcd-coefs
  807.                        p (math-poly-deriv-coefs p))))
  808.            (length p)))
  809.        (setq p (car temp))
  810.        (while (cdr temp)
  811.          (setq temp (nthcdr (nth 2 math-poly-modulus) temp)
  812.            p (cons (car temp) p)))
  813.        (and (setq temp (math-factor-poly-coefs p))
  814.         (math-pow temp (nth 2 math-poly-modulus))))
  815.       (t
  816.        (math-reject-arg nil "*Modulo factorization not yet implemented"))))
  817. )
  818.  
  819. (defun math-poly-deriv-coefs (p)
  820.   (let ((n 1)
  821.     (dp nil))
  822.     (while (setq p (cdr p))
  823.       (setq dp (cons (math-mul (car p) n) dp)
  824.         n (1+ n)))
  825.     (nreverse dp))
  826. )
  827.  
  828. (defun math-factor-contains (x a)
  829.   (if (equal x a)
  830.       1
  831.     (if (memq (car-safe x) '(+ - * / neg))
  832.     (let ((sum 0))
  833.       (while (setq x (cdr x))
  834.         (setq sum (+ sum (math-factor-contains (car x) a))))
  835.       sum)
  836.       (if (and (eq (car-safe x) '^)
  837.            (natnump (nth 2 x)))
  838.       (* (math-factor-contains (nth 1 x) a) (nth 2 x))
  839.     0)))
  840. )
  841.  
  842.  
  843.  
  844.  
  845.  
  846. ;;; Merge all quotients and expand/simplify the numerator
  847. (defun calcFunc-nrat (expr)
  848.   (if (math-any-floats expr)
  849.       (setq expr (calcFunc-pfrac expr)))
  850.   (if (or (math-vectorp expr)
  851.       (assq (car-safe expr) calc-tweak-eqn-table))
  852.       (cons (car expr) (mapcar 'calcFunc-nrat (cdr expr)))
  853.     (let* ((calc-prefer-frac t)
  854.        (res (math-to-ratpoly expr))
  855.        (num (math-simplify (math-sort-terms (calcFunc-expand (car res)))))
  856.        (den (math-simplify (math-sort-terms (calcFunc-expand (cdr res)))))
  857.        (g (math-poly-gcd num den)))
  858.       (or (eq g 1)
  859.       (let ((num2 (math-poly-div num g))
  860.         (den2 (math-poly-div den g)))
  861.         (and (eq (cdr num2) 0) (eq (cdr den2) 0)
  862.          (setq num (car num2) den (car den2)))))
  863.       (math-simplify (math-div num den))))
  864. )
  865.  
  866. ;;; Returns expressions (num . denom).
  867. (defun math-to-ratpoly (expr)
  868.   (let ((res (math-to-ratpoly-rec expr)))
  869.     (cons (math-simplify (car res)) (math-simplify (cdr res))))
  870. )
  871.  
  872. (defun math-to-ratpoly-rec (expr)
  873.   (cond ((Math-primp expr)
  874.      (cons expr 1))
  875.     ((memq (car expr) '(+ -))
  876.      (let ((r1 (math-to-ratpoly-rec (nth 1 expr)))
  877.            (r2 (math-to-ratpoly-rec (nth 2 expr))))
  878.        (if (equal (cdr r1) (cdr r2))
  879.            (cons (list (car expr) (car r1) (car r2)) (cdr r1))
  880.          (if (eq (cdr r1) 1)
  881.          (cons (list (car expr)
  882.                  (math-mul (car r1) (cdr r2))
  883.                  (car r2))
  884.                (cdr r2))
  885.            (if (eq (cdr r2) 1)
  886.            (cons (list (car expr)
  887.                    (car r1)
  888.                    (math-mul (car r2) (cdr r1)))
  889.              (cdr r1))
  890.          (let ((g (math-poly-gcd (cdr r1) (cdr r2))))
  891.            (let ((d1 (and (not (eq g 1)) (math-poly-div (cdr r1) g)))
  892.              (d2 (and (not (eq g 1)) (math-poly-div
  893.                           (math-mul (car r1) (cdr r2))
  894.                           g))))
  895.              (if (and (eq (cdr d1) 0) (eq (cdr d2) 0))
  896.              (cons (list (car expr) (car d2)
  897.                      (math-mul (car r2) (car d1)))
  898.                    (math-mul (car d1) (cdr r2)))
  899.                (cons (list (car expr)
  900.                    (math-mul (car r1) (cdr r2))
  901.                    (math-mul (car r2) (cdr r1)))
  902.                  (math-mul (cdr r1) (cdr r2)))))))))))
  903.     ((eq (car expr) '*)
  904.      (let* ((r1 (math-to-ratpoly-rec (nth 1 expr)))
  905.         (r2 (math-to-ratpoly-rec (nth 2 expr)))
  906.         (g (math-mul (math-poly-gcd (car r1) (cdr r2))
  907.                  (math-poly-gcd (cdr r1) (car r2)))))
  908.        (if (eq g 1)
  909.            (cons (math-mul (car r1) (car r2))
  910.              (math-mul (cdr r1) (cdr r2)))
  911.          (cons (math-poly-div-exact (math-mul (car r1) (car r2)) g)
  912.            (math-poly-div-exact (math-mul (cdr r1) (cdr r2)) g)))))
  913.     ((eq (car expr) '/)
  914.      (let* ((r1 (math-to-ratpoly-rec (nth 1 expr)))
  915.         (r2 (math-to-ratpoly-rec (nth 2 expr))))
  916.        (if (and (eq (cdr r1) 1) (eq (cdr r2) 1))
  917.            (cons (car r1) (car r2))
  918.          (let ((g (math-mul (math-poly-gcd (car r1) (car r2))
  919.                 (math-poly-gcd (cdr r1) (cdr r2)))))
  920.            (if (eq g 1)
  921.            (cons (math-mul (car r1) (cdr r2))
  922.              (math-mul (cdr r1) (car r2)))
  923.          (cons (math-poly-div-exact (math-mul (car r1) (cdr r2)) g)
  924.                (math-poly-div-exact (math-mul (cdr r1) (car r2))
  925.                         g)))))))
  926.     ((and (eq (car expr) '^) (integerp (nth 2 expr)))
  927.      (let ((r1 (math-to-ratpoly-rec (nth 1 expr))))
  928.        (if (> (nth 2 expr) 0)
  929.            (cons (math-pow (car r1) (nth 2 expr))
  930.              (math-pow (cdr r1) (nth 2 expr)))
  931.          (cons (math-pow (cdr r1) (- (nth 2 expr)))
  932.            (math-pow (car r1) (- (nth 2 expr)))))))
  933.     ((eq (car expr) 'neg)
  934.      (let ((r1 (math-to-ratpoly-rec (nth 1 expr))))
  935.        (cons (math-neg (car r1)) (cdr r1))))
  936.     (t (cons expr 1)))
  937. )
  938.  
  939.  
  940. (defun math-ratpoly-p (expr &optional var)
  941.   (cond ((equal expr var) 1)
  942.     ((Math-primp expr) 0)
  943.     ((memq (car expr) '(+ -))
  944.      (let ((p1 (math-ratpoly-p (nth 1 expr) var))
  945.            p2)
  946.        (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
  947.         (max p1 p2))))
  948.     ((eq (car expr) '*)
  949.      (let ((p1 (math-ratpoly-p (nth 1 expr) var))
  950.            p2)
  951.        (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
  952.         (+ p1 p2))))
  953.     ((eq (car expr) 'neg)
  954.      (math-ratpoly-p (nth 1 expr) var))
  955.     ((eq (car expr) '/)
  956.      (let ((p1 (math-ratpoly-p (nth 1 expr) var))
  957.            p2)
  958.        (and p1 (setq p2 (math-ratpoly-p (nth 2 expr) var))
  959.         (- p1 p2))))
  960.     ((and (eq (car expr) '^)
  961.           (integerp (nth 2 expr)))
  962.      (let ((p1 (math-ratpoly-p (nth 1 expr) var)))
  963.        (and p1 (* p1 (nth 2 expr)))))
  964.     ((not var) 1)
  965.     ((math-poly-depends expr var) nil)
  966.     (t 0))
  967. )
  968.  
  969.  
  970. (defun calcFunc-apart (expr &optional var)
  971.   (cond ((Math-primp expr) expr)
  972.     ((eq (car expr) '+)
  973.      (math-add (calcFunc-apart (nth 1 expr) var)
  974.            (calcFunc-apart (nth 2 expr) var)))
  975.     ((eq (car expr) '-)
  976.      (math-sub (calcFunc-apart (nth 1 expr) var)
  977.            (calcFunc-apart (nth 2 expr) var)))
  978.     ((not (math-ratpoly-p expr var))
  979.      (math-reject-arg expr "Expected a rational function"))
  980.     (t
  981.      (let* ((calc-prefer-frac t)
  982.         (rat (math-to-ratpoly expr))
  983.         (num (car rat))
  984.         (den (cdr rat))
  985.         (qr (math-poly-div num den))
  986.         (q (car qr))
  987.         (r (cdr qr)))
  988.        (or var
  989.            (setq var (math-polynomial-base den)))
  990.        (math-add q (or (and var
  991.                 (math-expr-contains den var)
  992.                 (math-partial-fractions r den var))
  993.                (math-div r den))))))
  994. )
  995.  
  996.  
  997. (defun math-padded-polynomial (expr var deg)
  998.   (let ((p (math-is-polynomial expr var deg)))
  999.     (append p (make-list (- deg (length p)) 0)))
  1000. )
  1001.  
  1002. (defun math-partial-fractions (r den var)
  1003.   (let* ((fden (calcFunc-factors den var))
  1004.      (tdeg (math-polynomial-p den var))
  1005.      (fp fden)
  1006.      (dlist nil)
  1007.      (eqns 0)
  1008.      (lz nil)
  1009.      (tz (make-list (1- tdeg) 0))
  1010.      (calc-matrix-mode 'scalar))
  1011.     (and (not (and (= (length fden) 2) (eq (nth 2 (nth 1 fden)) 1)))
  1012.      (progn
  1013.        (while (setq fp (cdr fp))
  1014.          (let ((rpt (nth 2 (car fp)))
  1015.            (deg (math-polynomial-p (nth 1 (car fp)) var))
  1016.            dnum dvar deg2)
  1017.            (while (> rpt 0)
  1018.          (setq deg2 deg
  1019.                dnum 0)
  1020.          (while (> deg2 0)
  1021.            (setq dvar (append '(vec) lz '(1) tz)
  1022.              lz (cons 0 lz)
  1023.              tz (cdr tz)
  1024.              deg2 (1- deg2)
  1025.              dnum (math-add dnum (math-mul dvar
  1026.                                (math-pow var deg2)))
  1027.              dlist (cons (and (= deg2 (1- deg))
  1028.                       (math-pow (nth 1 (car fp)) rpt))
  1029.                      dlist)))
  1030.          (let ((fpp fden)
  1031.                (mult 1))
  1032.            (while (setq fpp (cdr fpp))
  1033.              (or (eq fpp fp)
  1034.              (setq mult (math-mul mult
  1035.                           (math-pow (nth 1 (car fpp))
  1036.                             (nth 2 (car fpp)))))))
  1037.            (setq dnum (math-mul dnum mult)))
  1038.          (setq eqns (math-add eqns (math-mul dnum
  1039.                              (math-pow
  1040.                               (nth 1 (car fp))
  1041.                               (- (nth 2 (car fp))
  1042.                              rpt))))
  1043.                rpt (1- rpt)))))
  1044.        (setq eqns (math-div (cons 'vec (math-padded-polynomial r var tdeg))
  1045.                 (math-transpose
  1046.                  (cons 'vec
  1047.                        (mapcar
  1048.                     (function
  1049.                      (lambda (x)
  1050.                        (cons 'vec (math-padded-polynomial
  1051.                                x var tdeg))))
  1052.                     (cdr eqns))))))
  1053.        (and (math-vectorp eqns)
  1054.         (let ((res 0)
  1055.               (num nil))
  1056.           (setq eqns (nreverse eqns))
  1057.           (while eqns
  1058.             (setq num (cons (car eqns) num)
  1059.               eqns (cdr eqns))
  1060.             (if (car dlist)
  1061.             (setq num (math-build-polynomial-expr
  1062.                    (nreverse num) var)
  1063.                   res (math-add res (math-div num (car dlist)))
  1064.                   num nil))
  1065.             (setq dlist (cdr dlist)))
  1066.           (math-normalize res))))))
  1067. )
  1068.  
  1069.  
  1070.  
  1071. (defun math-expand-term (expr)
  1072.   (cond ((and (eq (car-safe expr) '*)
  1073.           (memq (car-safe (nth 1 expr)) '(+ -)))
  1074.      (math-add-or-sub (list '* (nth 1 (nth 1 expr)) (nth 2 expr))
  1075.               (list '* (nth 2 (nth 1 expr)) (nth 2 expr))
  1076.               nil (eq (car (nth 1 expr)) '-)))
  1077.     ((and (eq (car-safe expr) '*)
  1078.           (memq (car-safe (nth 2 expr)) '(+ -)))
  1079.      (math-add-or-sub (list '* (nth 1 expr) (nth 1 (nth 2 expr)))
  1080.               (list '* (nth 1 expr) (nth 2 (nth 2 expr)))
  1081.               nil (eq (car (nth 2 expr)) '-)))
  1082.     ((and (eq (car-safe expr) '/)
  1083.           (memq (car-safe (nth 1 expr)) '(+ -)))
  1084.      (math-add-or-sub (list '/ (nth 1 (nth 1 expr)) (nth 2 expr))
  1085.               (list '/ (nth 2 (nth 1 expr)) (nth 2 expr))
  1086.               nil (eq (car (nth 1 expr)) '-)))
  1087.     ((and (eq (car-safe expr) '^)
  1088.           (memq (car-safe (nth 1 expr)) '(+ -))
  1089.           (integerp (nth 2 expr))
  1090.           (if (> (nth 2 expr) 0)
  1091.           (or (and (or (> mmt-many 500000) (< mmt-many -500000))
  1092.                (math-expand-power (nth 1 expr) (nth 2 expr)
  1093.                           nil t))
  1094.               (list '*
  1095.                 (nth 1 expr)
  1096.                 (list '^ (nth 1 expr) (1- (nth 2 expr)))))
  1097.         (if (< (nth 2 expr) 0)
  1098.             (list '/ 1 (list '^ (nth 1 expr) (- (nth 2 expr))))))))
  1099.     (t expr))
  1100. )
  1101.  
  1102. (defun calcFunc-expand (expr &optional many)
  1103.   (math-normalize (math-map-tree 'math-expand-term expr many))
  1104. )
  1105.  
  1106. (defun math-expand-power (x n &optional var else-nil)
  1107.   (or (and (natnump n)
  1108.        (memq (car-safe x) '(+ -))
  1109.        (let ((terms nil)
  1110.          (cterms nil))
  1111.          (while (memq (car-safe x) '(+ -))
  1112.            (setq terms (cons (if (eq (car x) '-)
  1113.                      (math-neg (nth 2 x))
  1114.                    (nth 2 x))
  1115.                  terms)
  1116.              x (nth 1 x)))
  1117.          (setq terms (cons x terms))
  1118.          (if var
  1119.          (let ((p terms))
  1120.            (while p
  1121.              (or (math-expr-contains (car p) var)
  1122.              (setq terms (delq (car p) terms)
  1123.                    cterms (cons (car p) cterms)))
  1124.              (setq p (cdr p)))
  1125.            (if cterms
  1126.                (setq terms (cons (apply 'calcFunc-add cterms)
  1127.                      terms)))))
  1128.          (if (= (length terms) 2)
  1129.          (let ((i 0)
  1130.                (accum 0))
  1131.            (while (<= i n)
  1132.              (setq accum (list '+ accum
  1133.                        (list '* (calcFunc-choose n i)
  1134.                          (list '*
  1135.                            (list '^ (nth 1 terms) i)
  1136.                            (list '^ (car terms)
  1137.                              (- n i)))))
  1138.                i (1+ i)))
  1139.            accum)
  1140.            (if (= n 2)
  1141.            (let ((accum 0)
  1142.              (p1 terms)
  1143.              p2)
  1144.              (while p1
  1145.                (setq accum (list '+ accum
  1146.                      (list '^ (car p1) 2))
  1147.                  p2 p1)
  1148.                (while (setq p2 (cdr p2))
  1149.              (setq accum (list '+ accum
  1150.                        (list '* 2 (list '*
  1151.                                 (car p1)
  1152.                                 (car p2))))))
  1153.                (setq p1 (cdr p1)))
  1154.              accum)
  1155.          (if (= n 3)
  1156.              (let ((accum 0)
  1157.                (p1 terms)
  1158.                p2 p3)
  1159.                (while p1
  1160.              (setq accum (list '+ accum (list '^ (car p1) 3))
  1161.                    p2 p1)
  1162.              (while (setq p2 (cdr p2))
  1163.                (setq accum (list '+
  1164.                          (list '+
  1165.                            accum
  1166.                            (list '* 3
  1167.                              (list
  1168.                               '*
  1169.                               (list '^ (car p1) 2)
  1170.                               (car p2))))
  1171.                          (list '* 3
  1172.                            (list
  1173.                             '* (car p1)
  1174.                             (list '^ (car p2) 2))))
  1175.                  p3 p2)
  1176.                (while (setq p3 (cdr p3))
  1177.                  (setq accum (list '+ accum
  1178.                            (list '* 6
  1179.                              (list '*
  1180.                                (car p1)
  1181.                                (list
  1182.                                 '* (car p2)
  1183.                                 (car p3))))))))
  1184.              (setq p1 (cdr p1)))
  1185.                accum))))))
  1186.       (and (not else-nil)
  1187.        (list '^ x n)))
  1188. )
  1189.  
  1190. (defun calcFunc-expandpow (x n)
  1191.   (math-normalize (math-expand-power x n))
  1192. )
  1193.  
  1194.  
  1195.  
  1196.